Metabolic model-based ecological modeling for probiotic design

The microbial community composition in the human gut has a profound effect on human health. This observation has lead to extensive use of microbiome therapies, including over-the-counter ‘probiotic’ treatments intended to alter the composition of the microbiome. Despite so much promise and commercial interest, the factors that contribute to the success or failure of microbiome-targeted treatments remain unclear. We investigate the biotic interactions that lead to successful engraftment of a novel bacterial strain introduced to the microbiome as in probiotic treatments. We use pairwise genome-scale metabolic modeling with a generalized resource allocation constraint to build a network of interactions between taxa that appear in an experimental engraftment study. We create induced sub-graphs using the taxa present in individual samples and assess the likelihood of invader engraftment based on network structure. To do so, we use a generalized Lotka-Volterra model, which we show has strong ability to predict if a particular invader or probiotic will successfully engraft into an individual’s microbiome. Furthermore, we show that the mechanistic nature of the model is useful for revealing which microbe-microbe interactions potentially drive engraftment.


Introduction
Microbiome research has come to encompass key areas of disease, ranging from infections (Antharam et al., 2013;Honda and Littman, 2012;Battaglioli et al., 2018) and cancer prevention (Moss and Blaser, 2005;Walther-António et al., 2016;Kim et al., 2020) to systemic immune and neurological responses (Severance et al., 2016;Kang et al., 2014;Chen et al., 2016).The effect of the microbiome on health is now undeniable, and every year in the US over 400,000 people collectively spend $1 billion dollars on over-the-counter probiotics intended to alter their microbiome (Kristensen et al., 2016).Many of the purported interactions between microbes and health involve resident microbiota and their interactions with the host, i.e., the interface between microbial ecology and human health.The goal of microbiome-targeted interventions is therefore to promote health by "restoring and maintaining the microbiota and the crucial health-associated ecosystem services that it provides" (Costello et al., 2012).
Despite the many links between the microbiome and health, our ability to deploy probiotics to modify the microbiome as intended has been met with relatively little success (Mullard, 2016;Zhu et al., 2019;Yuan et al., 2017;Zhao et al., 2021;Wu et al., 2017).Studies looking at the ecological effects of probiotic administration show that administration of a probiotic is not sufficient to alter the community in the desired way.Specifically, engraftment of the administered microbial species is often limited, with only one-third to one-half of patients showing any signs of medium or long-term engraftment (Maldonado-Gómez et al., 2016;Pudgar et al., 2021).We offer the argument that probiotic interventions are primarily ecologic in nature; their purpose is to reshape the complex microbial communities in our body in beneficial ways.Therefore, to predict whether a probiotic has the desired effects in the gut microbial community, we need more studies examining the ecology of probiotic interventions (Walter et al., 2018).A mechanistic, personalized approach to probiotic design-one rooted in empirical metabolic data and ecologic principles-has the potential to propel the field forward.
Previous work trying to predict engraftment has been mostly restricted to fecal microbiome transplant (FMT) studies using non-mechanistic classifiers (Smillie et al., 2018;Podlesny et al., 2021).While potentially predictive, such classifier approaches are sensitive to the underlying assumptions or conditions in which the study is carried out.Because such models are built using statistical methods on data that is assumed to be uniformly collected, these models cannot be generalized to new circumstances.It would be unexpected to see predictions built from patients with diarrhea, undergoing bowel prep, taking antibiotics, and given an FMT accurately predict what happens in patients taking orally administered probiotics.
The goal of this work is to examine the use of metabolic modeling-informed population dynamic approaches.This builds on top of related work in both constraint-based metabolic modeling and population models such as Lotka-Volterra.It is worth highlighting that the use of dynamic flux balance analysis for population models has also been a well-published approach.Despite these successes, there are a number of practical drawbacks for communities of high complexity such as labor-intensive interpretation and high computational complexity.
Population models such as the generalized Lotka-Volterra model are popular tools for understanding microbial community dynamics in a mechanistic manner (Stein et al., 2013;Friedman et al., 2017;Angulo et al., 2019;Kuntal et al., 2019).However, these models are in general difficult to parameterize, with state of the art gradient-matching procedures requiring somewhat dense time-longitudinal data with many replicates Bucci et al. (2016).Furthermore, the authors have previously shown that parameters fit from data to these models do not extend to novel environmental situations, and may even change with the addition of a new taxa to the community (Brunner and Chia, 2019).These drawbacks make such mechanistic population models impractical for predicting engraftment.However, by leveraging metabolic modeling we are able to parameterize population models in a way that can be easily adapted to novel environments and does not require dense time-longitudinal data.This allows us to use these models to predict microbial engraftment into a community.See fig. 1 for a comparison of our method with standard parameter fitting.
In this paper, we present a method to predict engraftment of an invader into a microbial community in the following manner.First, we construct an interaction network of the microbial taxa found in a sample from the community using pairwise flux balance analysis with resource allocation constraints (Kim et al., 2022) with genome-scale metabolic models included in the AGORA database (Magnúsdóttir et al., 2017).We then make a prediction based off of one or more of six dynamical systems models parameterized by this network-the generalized Lotka-Volterra (LV) model Edelstein- Keshet (2005); Stein et al. (2013); Friedman et al. (2017), an adjustment to the LV model we call the "antagonistic Lotka-Volterra" model (AntLV), another adjustment to the LV model we call the "inhibitory Lotka-Volterra" model (InhibtLV), a fourth non-linear model based on the replicator equation Madec and Gjini (2020); Gjini and Madec (2021) which is similar to the LV model, and two similar linear models based on node balancing (NodeBalance) and random walks (Stochastic) (see fig. 2).
We test the predictive potential of each dynamical model by predicting the outcomes of microbiome invasion experiments (Battaglioli et al., 2018;Maldonado-Gómez et al., 2016;Huang et al., 2021) from the initial presence/absence of species in each sample.We choose to test all six models because no definitive "best model" has been established.The LV model is widely used, but we found that this model could lead to uncontrolled simulated growth.We developed the AntLV and Figure 1.Population models such as Lotka-Volterra generally require dense time-longitudinal data to accurately parameterize.In this work, we leverage genome-scale metabolic modeling to parameterize population models with only genomic data from a single time-point.We are unable to parameterize these models using standard techniques due to the sparsity of the data.
InhibitLV models in order dampen the numerical instability of the LV model.We also include the replicator model because has recently been shown to provide insight into microbial population dynamics (Madec and Gjini, 2020).Finally, we include the two linear models because they provide a substantial reduction in computational complexity, and so are more practical to use as long as they provide adequate predictive power.
Because the specific invasion experiments involved invader strains not present in the AGORA database, we repeated the experiment for each species-level match to the invader in the database.We show that this method has good predictive potential depending on the choice of dynamical system used to score the sub-graphs and choice of AGORA database match to the invader strain.Furthermore, we perform two types of sensitivity analysis to demonstrate that the mechanistic nature of the model provides additional insight into the impact of the various components of the network.We perform simulated knock-out experiments to test how sensitive engraftment is to each community member, and we use parameter sensitivity analysis to test how sensitive engraftment is to each network connection.

Results
We first examine the ability of our metabolic modeling-based approach to successfully predict engraftment versus non-engraftment for different microbial species introduced orally across different experimental or clinical trial settings.The three studies we use to test predictive power all involved the introduction of an invading taxa into an established microbial community.Two of these studies examine the introduction of candidate probiotics, B. longum and L. plantarum, respectively, into human subjects (Maldonado-Gómez et al., 2016;Huang et al., 2021).The third examines the introduction of the pathogen C. difficile (Battaglioli et al., 2018) into mice that had been inoculated with microbial communities from healthy and disbiotic human donors.Note that while this study introduces a pathogen rather than a probiotic, the ecological principle of engraftment into an existing community is the same.All three studies include microbiome samples from prior to introduction of the invader, which we use to generate predictions, as well as samples from Predictions are made for each sample using 6 different dynamical models Four quadratic models.after the community was allowed to reform, which we use to score our predictions.
As a metric of classification success for the C. difficile and B. longum data-sets, we use the area under the curve of the receiver-operator characteristic (AUC-ROC).This metric provides a measure of performance based on the model's ability to identify true positives while avoiding false positives, so that 1 is perfect classifier performance and 0.5 is equivalent to random classification (i.e.flipping a coin for each sample).We use Kendall-tau rank correlation as a metric of classifier success for the L. plantarum data-set, comparing our test-score with the observed abundance.This is necessary because AUC-ROC requires binary classification of the test data, which was unavailable for this data-set.
In figs. 3 to 5, we report the AUC-ROC of each dynamical model for each experiment, as well as the AUC-ROC for support vector machine and random forest classification or regression (see fig. 2 for the dynamics used).Tables of these results and associated estimated p-values can be found in the supplementary file supp_tables.pdf,(results in Tables 1,3 and 6, p-values in Tables  2,4,7, and standard machine learning's performance in Tables 5 and 8).Study data from (Battaglioli et al., 2018) displayed clear differences between engrafter and non-engrafter samples in both species composition and -diversity (see Supplementary Figure 3 in supp_figures.pdf).For this reason, any reasonable classifier should be expected to differentiate the two groups.Indeed, our method produced perfect classification on the C. difficile data-set for every of choice of dynamics or GEM representing the invading C. difficile.We may consider this performance a necessary, but not sufficient criteria for the quality of the method.
On the other hand, data from the studies (Maldonado-Gómez et al., 2016;Huang et al., 2021) was much more muddled.In both of these studies, our method generally outperforms the classical machine-learning techniques to which we compared it.The Lotka-Volterra dynamics performed the best, especially when modified to be made antagonistic or inhibitory.Inspection of individual simulations revealed that this improvement is due to those modifications preventing finite-time blow up of solutions.
We tested three general types of community model that can be associated with a network of interactions: generalized Lotka-Volterra dynamics (gLV, antagonistic LV, and inhibitory LV), replicator dynamics, and linear dynamics (node balancing and stochastic).The two linear models we used were both based on balance within a network; one based on balancing flux through the network nodes and the other based on the equilibrium of a stochastic system.Neither model performed well, suggesting that the low complexity of linear modeling is not sufficient for prediction.We observe much better predictive power from the two types of quadratic models used.The main difference between the generalized Lotka-Volterra model and the replicator model is that the replicator model includes a "community effect" that is equal for every member of the community.Although this effect has been shown to be useful in predicting community composition (Gjini and Madec, 2021;Madec and Gjini, 2020), it does not improve predictive power of our method on the data that we tested.

Interpreting Mechanism -Sensitivity to Hyper-parameters
One major advantage to any mechanistic model is that we may measure the sensitivity of our results to perturbations in the model.Here, we investigate the results of perturbing the model of the B. longum experiments in two ways -simulating "knock-out" experiments and computing sensitivity to changes in interaction strength.
First, we simulate knock-out experiments by removing a taxa from the full AGORA network.As a result, this organism will be removed from any sample it was previously present in.We then measure the effect this has on our prediction of B. longum's likelihood of engrafting based on the antagonistic Lotka-Volterra dynamics.
Next, we measure the sensitivity of B. longum growth to each interaction parameter in the antagonistic Lotka-Volterra model.That is, we measure the effect of perturbing each parameter individually on the simulated abundance of B. longum at equilibrium according to antagonistic Lotka-Volterra dynamics.

Sensitivity to Community Members
In order to investigate how sensitive B. longum engraftment is to each of the other species present in any sample of the B. longum experimental data set, we simulate "knock-out" experiments and observe the effect this has on our prediction of engraftment.We used the model of B. longum infantis 157F NC to model invading B. longum because this choice provided the best predictive value when using the full AGORA network and antagonistic Lotka-Volterra dynamics (see Table 3 of supp_tables.tex).For each simulated knock-out, we removed a taxa from every sample it was present in and repeated our predictive procedure, using only antagonistic Lotka-Volterra dynamics for our prediction.We recorded the difference in sample score for engraftment of B. longum (eq.( 5)) as well as the change in AUC-ROC of our set of predictions.
We present the five knock-outs that had the largest average effect on our sample score across all samples of the B. longum data set in table 1.It is interesting to note the difference that each knock-out has on the predictive power of the model.For instance, removing L. delbrueckii subsp.bulgaricus ATCC BAA 365, which on average increases the likelihood of predicting that B. longum will engraft, causes our model to have very poor predictive value (AUC-ROC of 0.40).This suggests that detection of L. delbrueckii is important in prediction.Interestingly, this effect is due to a difference in knock-out effects between actual engrafters and actual non-engrafters.On average, L. delbrueckii subsp bulgaricus ATCC BAA 365 knock-out increased the engraftment score of B. longum by 0.086 in non-engrafting samples and only 0.02 in engrafting samples.
This experiment demonstrates that incorporating mechanism into prediction can provide useful insight beyond prediction.Because our model considers the network of interactions between microbial taxa, we are able to experiment with the effects of each individual taxa by using simulated experiments like knock-outs.The "Sample Proportion" column gives the proportion of samples in the data-set that contain the organism that was knocked out.The final column shows the impact on our predictions of removing the organism from the analysis.

Sensitivity to Interactions
Our method is based on the network of interactions between microbial species, which ties together the individual interactions between each pair of species into one complex picture.This means that the interaction between two species unrelated to B. longum may have an effect B. longum's growth.
We can compute this effect by computing how B. longum's simulated abundance changes as we vary each parameter.Precisely, we can compute the derivative of B. longum with respect to each parameter in the model.We once again limit this analysis to B. longum infantis 157F NC.
We compute the sensitivity of B. longum's growth according to antagonistic Lotka-Volterra dynamics to each interaction parameter directly using the chain rule (see, for example (Zi, 2011)).This has the following form: (1) Equation ( 1) allows us to solve a system of differential equations for each parameter in the network to determine the sensitivity of invader growth to that parameter.We present the results of this analysis, averaged across each sample in the data set, in table 2. It is unsurprising that the growth of B. longum was most sensitive to two direct regulations in the network, i.e. regulations that have B. longum as their target.However, we observe that B. longum is also sensitive to the regulation of G. pamelaeae by B. crossotus, as well as a number of other regulations that are mediated by B. crossotus.These indirect effects suggests that B. crossotus plays a large role in B. longum growth even if it does not directly regulate B. longum.
This analysis demonstrates one advantage to network based prediction.Networks synthesize sets of individual interactions into a cohesive whole, and in doing so reveal indirect interactions that may be crucial to predicting microbial community composition.

Discussion
Our work exams the scenario where the goal of probiotic intervention is a long-lasting alteration of the resident host microbiota.Given that the links between microbes and health have a basis in observations of the resident human microbiota (Gupta et al., 2020), one common assumption is that alterations to that resident microbiota would impact health outcomes.Along those lines, we demonstrate that genome-scale metabolic modeling can be used in a simple way to predict the outcome of species invasion experiments with minimal study data required.This suggests that genome-scale models (GEMs) can provide value in understanding microbial community dynamics from cross-sectional microbiome relative abundance profiling.
Genome-scale metabolic modeling uses genomic data to predict the growth and resource use of a microbial population by considering the entire internal metabolism of the species.nique can be extended to community modeling in a number of ways (Zomorrodi and Maranas, 2012;Chan et al., 2017;Diener et al., 2020;Kim et al., 2022), all with relatively high levels of complexity.Community dynamics can be inferred from GEMs using a model known as dynamic flux balance analysis (FBA), which uses a GEM for each taxa in the community to infer growth rates and dynamic resource usage (Brunner and Chia, 2020).While dynamic FBA provides a complete picture of community dynamics according to GEMs, it represents a very complex model that is difficult to analyze and simulate, and still contains a set of unknown parameters.Here, we simplify the modeling considerably by assuming that these underlying dynamics lead to a network of emergent interactions between microbes, which we can infer from pairwise genome-scale modeling.
The generalized Lotka-Volterra model and other similar microbial network models are popular tools for understanding community dynamics (Friedman et al., 2017;Fisher and Mehta, 2014;Angulo et al., 2019).Such models are difficult to parameterize for a myriad of reasons.Good parameterization requires relatively dense time-longitudinal data of absolute, rather than relative, community abundance (Bucci et al., 2016;Kuntal et al., 2019).By taking advantage of genomescale metabolic models (GEMs), we provide a parameterized network that does not require any time-longitudinal data.In fact, our method as presented here makes use of previously published GEMs, meaning that only binary presence/absence information is necessary for a prediction.Additionally, GEMs can be built directly from the genomes found in a sample using automated tools such as CarveMe (Machado et al., 2018), meaning that our method can be extended to include taxa not found in the AGORA database.
Another strength of our approach is it's interpretability.The B. longum engraftment analysis provides us with an example motivation, which would be to identify other partner microbes that could improve the responsiveness to the engraftment of a target probiotic.In this instance, there is a very strong link between L. delbrueckii and B. longum engraftment in our model.This is consistent with co-culture studies that show L. delbrueckii and B. longum tolerate stressful conditions better together than when they are alone (Sánchez et al., 2013), suggesting some type of cooperation between the species.While we cannot know the mechanism of this cooperation from literature, our modeling suggests that it is at least reflected in part by metabolism.Similarly intriguing is the cooperative link between B. crossotus and B. longum given the known cross-feeding interactions between many other butyrate-producing organisms (Falony et al., 2006;De Vuyst and Leroy, 2011;Rios-Covian et al., 2015)-a metabolite that Butyrivibrio produce in large amounts (Cotta and Forster, 2006).
Lastly, we note that we have previously shown that species-species interaction modeling, i.e. models built from interactions between microbes, do not capture the complexities of microbial community dynamics that emerge as communities change in composition (Brunner and Chia, 2019).
Here, we mitigate this shortcoming by determining interaction parameters from pairwise models under specific metabolic conditions, providing limited environmental context to our method.However, it is unlikely that these interactions remain constant as the microbial community manipulates its environment.We conjecture that prediction can be improved by accounting for changes in microbial interaction as the environment changes.In upcoming work, we demonstrate that dynamic FBA implies that a microbial community behaves according to a discrete sequence of interaction networks over time.Incorporating this dynamic behavior may improve prediction with only a modest increase in model complexity, as long as this sequence of networks can be efficiently determined.

Inferring interaction parameters
To make predictions without fitting any model to time series data, we use a network of interactions implied by a technique known as flux balance analysis (FBA).FBA is a technique used to predict growth rates of microbes using genome-scale information about their internal metabolisms (Lewis et al., 2012).This technique requires a genome-scale metabolic model (GEM) which represents the set of metabolic pathways of a microbe, as well as information about the environmental metabolites available.An optimization problem is then solved to predict a growth rate of the microbe.Similarly, a technique known as "joint flux balance analysis" (Zomorrodi and Maranas, 2012;Chan et al., 2017;Kim et al., 2022) can be used to estimate growth rates of pairs of microbes using a GEM from each microbe.
We use flux balance analysis with a "western diet" environmental condition to simulate growth of the 818 taxa included in the AGORA database of GEMs (Magnúsdóttir et al., 2017).Next, we use joint flux balance analysis to simulate the growth of every pair of these 818 microbes.We use joint FBA with a set of resource allocation constraints, which have been shown to improve predictions in pairwise experiments (Kim et al., 2022).We construct our interaction network by taking a complete, weighted, directed graph with the 818 taxa as nodes and the log-ratio of simulated growth in pairs and alone as edge weights.That is, if simulated growth of species alone is , and growth of species when coupled in a pair with species is , we weight the edge from species to species as * = log The log-ratio is used for simplicity and ease of computation.We also re-scale our network so that all edge weights are in the interval [−1, 1], meaning we take as our edge weights This rescaling is done for computational convenience, and can be viewed as a time-rescaling of the network dynamics which does not effect our predictions.We also trim the graph of the weaker edges to reduce the effect of spurious relationships between metabolic models.We tested our method using the strongest % of edges in absolute value for ranging from 10 to 100.Predictions for C. difficile and L. platarum did not show any pattern of effect of , while predictions for B. longum improved or stayed level as was decreased to some threshold, roughly = 25, after which results became less consistent (see supplemental Figures 7-15 in supp_figures.pdf).Here, we present results using = 25, meaning we keep only the strongest quartile of edges in absolute value.The result is a network of 818 nodes and 167281 edges, with a mean absolute edge weight of 0.515 and variance in absolute edge weight of 0.0382.We will refer to this final network as the "AGORA network".

Creating induced sub-graphs
In order to carry out our analysis using a network constructed from the AGORA database of GEMs (Magnúsdóttir et al., 2017), it was necessary to match the taxa present in the data with GEMs in the database.The B. longum study originally published in (Maldonado-Gómez et al., 2016) made available metagenomic sequence data, while the other two studies made taxonomic profiles available.We created taxonomic profiles for the B. longum study using Bracken (Lu et al., 2017).
For each data-set, we matched taxonomic profiles to the set of AGORA models by querying the NCBI taxonomy database for nearest match taxa and comparing to the known taxonomy of the AGORA models.Only taxa with a species-level match to an AGORA model were included in our analysis.The result is that our analysis was performed using only a subset of the taxa in each sample.For the C. difficile study published in (Battaglioli et al., 2018), we used an average of 36.2% of the relative abundance of each sample; for the B. longum study published in (Maldonado-Gómez et al., 2016), this was an average of 8.42% of the relative abundance; and for the L. plantarum study published in (Huang et al., 2021), this was 14.8% of the relative abundance.

Assessing Receptivity
In order to assess how receptive a sub-graph is to an invading taxa, we use six different dynamical systems that can be associated with a pairwise graph, along with a composite score using all six.Each dynamical system allows us to simulate the community to equilibrium and score the performance of the invading taxa by simulated final abundance or time to extinction.The composite score can then be calculated simply as an average of the scores from the six dynamical systems.The dynamical systems fell into three categories -generalized Lotka-Volterra models, replicator equation-based models, and linear models.

Generalized Lotka-Volterra Models
The generalized Lotka-Volterra (gLV) model (Edelstein-Keshet, 2005) of a community of species is written as follows: Notice that we use a version of this model that assumes an intrinsic growth rate of 1 for each taxa because fitting accurate intrinsic growth rates would require additional data, while we wish to assess the performance of our methods without the need for parameter fitting.We use three variations of the generalized Lotka-Volterra models, based on adjustment to the edge weights.These are: (a) a model with the network edge weights as interaction parameters, = , which we refer to as "Lotka-Volterra dynamics", (b) a model with alll edge weights shifted to be negative, = ( −1) ∕2, which we refer to as "antagonistic Lotka-Volterra dynamics", and (c) a model with added self-inhibition, = for ≠ and = −1, which we refer to as "inhibitory Lotka-Volterra dynamics".Models (b) and (c) are adjustments to the generalized Lotka-Volterra model that we implemented to prevent uncontrolled growth in simulation.In simulation, the standard version (a) exhibited blowup of solutions which slowed computation and worsened prediction.Model (b) prevents blow-up by assuming that all interactions must be negative, or antagonistic, and adjusts all interactions so their relative effects are the same but absolute effects are negative.Model (c) prevents blow-up by using logistic growth for each individual microbe in isolation, rather than exponential growth in the standard gLV model.
The gLV model can display an array of behaviors, including multi-stability and chaos.For this reason, we use a Monte-Carlo sampling approach to generating predictions, taking repeated random draws of initial conditions and simulating forward.In our experiments, we use a number of draws equal to the number of nodes in the network.For each trial, we score the network based on the simulated relative abundance of the invader at some large time .It is possible that the invading taxa either dominates the community (i.e.relative abundance approaches 1) or becomes extinct in our simulation.We take this into account by setting a score as follows: = 1 3 time to extinction + relative abundance at time + time to community dominance (5) for each trial and averaging over all trials.

The Replicator Model
Recently, the replicator equation has been proposed as a model to study microbial community dynamics (Madec and Gjini, 2020;Gjini and Madec, 2021).This equation takes the form where is the matrix of .We take the network edges weights as interaction parameters, = and refer to this model as "replicator dynamics".As with the generalized Lotka-Volterra model, this model can display an array of complex behaviors, so we again use a Monte-Carlo sampling approach, scoring each sample according to eq. ( 5) and averaging over all samples.
Linear Models Finally, we use two linear models that arise from considering the network as a mathematical object rather than considering the biological reality of the microbial community.We include these models for two reasons.The first is that they are fundamental to the network structure of the community, and the second is that they are much more computationally efficient and so convenient to use.
The first of these models, which we refer to as "node balancing dynamics" is based on the idea of balancing a flow of mass into and out of each node.Precisely, we rescale the network edge weights to be in the interval [0, 1], taking ̂ = + 1 2 because the notion of node-balancing requires only positive interactions.Next, we compute the following graph Laplacian.Let be the matrix of rescaled edge weights ̂ so that species and are in the sample being scored.Then the Laplacian is defined as where the diagonal matrix is the weighted degree matrix of the re-scaled networks, = ∑ ̂ .We then score the network based on the dynamical model = . (8) Letting be the index in the vector corresponding to the invading taxa, we can compute a score for the network using the dominant eigenvector (i.e. the eigenvector corresponding to the largest eigenvalue) of .That is, if is the dominant eigenvector of , we take = as the score the network.
The second linear model is very closely related to the first, and is based on the idea of a Markov process jumping between the nodes of the network and so we refer to this model as "stochastic dynamics".Again, let be the matrix of re-scaled edge weights ̂ .From this, we construct a stochastic matrix ̂ by re-scaling so that ̂ = ∑ -the row sums of ̂ are 1.̂ now represents the transition matrix of a Markov jump process that moves between the nodes of the network, and we score the network according to the equilibrium probability that the process is in the node corresponding to the invading taxa.This can be calculated by computing the left eigenvectors of the matrix ̂ and noting the property of stochastic matrices that there exists a left eigenvalue = 1 and any corresponding eigenvector is a stationary distribution for the process.Precisely, if 1 is the left eigenvector corresponding to the eigenvalue = 1 and is the index of the invading taxa, we score the network = 1 .In brief, we generate an interaction network of 808 AGORA models using pairwise joint flux balance analysis.To produce a prediction of engraftment for a given sample, we use the taxa present in the sample to generate an induced sub-graph of the full AGORA network.This is then used in a dynamical system to generate a prediction of engraftment.

Prediction and evaluation
We  2016) separated samples into successful and unsuccessful invasions, allowing us to use our method as a classifier.For these studies, we were able to compute scores for each sample, and vary the discrimination thresh-hold for the binary prediction across the observed values.From this, we could compute a ROC curve and integrate (commonly referred to as the "area under the curve" of the ROC, or AUC-ROC).AUC-ROCs take values in the interval [0, 1] and, in general, an AUC-ROC greater than 0.5 indicates positive predictive value of the model.We compared the AUC-ROC across the choices of dynamics and further compared to classification using support vector machine (SVM) and random forest (RF) classification (Pedregosa et al., 2011).The SVM and RF classifications were performed using both the relative abundances from the data as well as binary (presence/absence) forms of the data, which matches our method.
The study Huang et al. (2021) did not classify samples, and so we compared our predictions to the final time-point of that study using two well known measures of rank correlation: Kendalltau rank correlation and Spearman rank correlation.Again, we compared across the choices of dynamics and to support vector machine regression and random forest regression (Pedregosa et al., 2011).
For all three studies, we estimated the significance of our predictions against predictions from a null-model created by replacing the AGORA network with a random network.To match the AGORA network in distribution, we used an approximate edge-swapping procedure.Precisely, we con-structed a complete graph on the 818 AGORA models by drawing edge weights from the AGORA network, with replacement.We then filtered the network as in the construction of the AGORA network by keeping only the top 25% of edges in absolute strength.We then computed AUC-ROC and rank correlations for predictions generated by this null model.We repeated this experiment 400 times in order to estimate a distribution of predictions from the null model.
Because we did not have strain specific models for the invading species in each data-set, we repeated our analysis for each model that matched at the species level with the invader.
exponentially modified by positive and negative interactions and a community wide interaction term determined by joint FBA Two linear models.Node Balance Biomass flows linearly through the network with conductance determined by joint FBA.Stochastic Biomass behaves as a random walk with transition probabilities determined by joint FBA.

Figure 2 .
Figure 2. We used six dynamical systems to assess the engraftment potential of an invading taxa.All six are parameterized based on a set of network connections derived from metabolic modeling.The four quadratic models are variations on the generalized Lotka-Volterra model.The two linear models are related to diffusion on the graph.
Figure 3. (A) All network dynamics (blue bars) tested for the four AGORA C. difficile models and the standard machine learning techniques (red bars) accurately predicted the engraftment versus non-engraftment as determined by the mean area under the receiver operator characteristic curve (AUC-ROC) (with 1 being perfect classification of engraftment and 0.5 being random chance).LV = Lotka-Volterra, AntLV = antagonistic Lotka-Volterra, inhibitLV = self-inhibition Lotka-Volterra, Replicator = replicator equation, NodeBalance = linear node balancing, Stochastic = random walk stochastic model, Composite = average of the previous six, SVM (Cont.)= support vector machine with continuous data, SVM (Disc.)= support vector machine with discretized data, RF (Cont.)= random forest with continuous data, RF (Disc.)= random forest with discretized data.(B) Example receiver operator characteristic curve (ROC) for a C. difficile genome visualizing the perfect prediction of engraftment -all true positives are identified before any false positives.(C) Example comparison between the antagonistic Lotka-Volterra model (antLV) predictions of C. difficile CD196 engraftment with the network parameters from joint FBA (red line) versus 400 null model networks (blue histogram).Null model networks consist of random interaction parameters generated using a label-swapping procedure to preserve the distribution of edge weights from the AGORA network.The clear superior performance of engraftment predicted with the AGORA network versus random networks indicates some level of biological significance.

Figure 6 .
Figure6.Schematic of the modeling process.In brief, we generate an interaction network of 808 AGORA models using pairwise joint flux balance analysis.To produce a prediction of engraftment for a given sample, we use the taxa present in the sample to generate an induced sub-graph of the full AGORA network.This is then used in a dynamical system to generate a prediction of engraftment.
used data from Battaglioli et al. (2018); Maldonado-Gómez et al. (2016); Huang et al. (2021) in order to test the predictive power of the method.Data from Battaglioli et al. (2018) consisted of microbial communities taken from previously germ-free mice that had been inoculated with microbial communities from one of two dysbiotic or one of two healthy donors.In that study, these mice were then exposed to Clostridioides difficile to confirm that the dysbiotic mice were more susceptible to engraftment.Data from Maldonado-Gómez et al. (2016) consisted of bacterial community composition of fecal samples taken over the course of experiments in which Bifidobacterium longum AH1206 was administered as an oral probiotic to 23 subjects.Subjects were then differentiated into "engrafters" and "non-engrafters" based on the survival of the probiotic strain.Finally, data from Huang et al. (2021) consisted of bacterial community composition from fecal samples from 12 healthy participants who were administered Lactiplanitbaccilus plantarum HNU082 as an oral probiotic.In this study, we differentiated invasion success using the final time-point of each experiment.Studies Battaglioli et al. (2018); Maldonado-Gómez et al. (

Table 1 .
Top 5 most impactful knock-out experiments by change in sample score (Avg.Score Difference).

Table 2 .
This tech-Sensitivity of B. longum growth to the 10 most impactful interaction parameters in the AGORA network (excluding self-loops).